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FINAL REPORT 

The objective of the work over the past year was to develop and use models of 
the lithium/thionyl chloride (Li/SOC^) cell to 

1. aid in understanding the physical phenomena occurring in the cell, 

2. reduce time intensive and costly experimental work, and 

3. expedite the finding of acceptable and safe, yet optimal, designs. 

These objectives have been partially satisfied by parameter estimation work with a 
one-dimensional mathematical model, developed earlier, and the heat flow patterns 
obtained from a two-dimensional thermal model. 

This final report describes the progress which has been made in modeling the 
lithium/thionyl chloride cell over the past year. The one-dimensional mathematical 
model for a lithium/thionyl chloride cell, developed earlier, was used in conjunction 
with a parameter estimation technique to show how useful cell parameters could 
be obtained from experimental data. A paper on this work has been submitted 
to the Journal of the Electrochemical Society for publication and can be found in 
appendix A. A simplified two-dimensional model for the spirally wound design of 
this battery was developed and used to investigate the heat flow within the cell. A 
paper describing this work has been accepted for publication in the Journal of the 
Electrochemical Society and is included in this report as appendix B. 
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ABSTRACT 


A one-dimensional mathematical model for the lithium/thionyl chloride primary 
cell is used in conjunction with a parameter estimation technique in order to estimate 
the electro-kinetic parameters of this electrochemical system. The electro-kinetic 
parameters include the anodic transfer coefficient and exchange current density of the 
lithium oxidation, a a ,i and the cathodic transfer coefficient and the effective 

exchange current density of the thionyl chloride reduction, o C)2 and a°t 0 ,2,r«/i and a 
morphology parameter, £. The parameter estimation is performed on simulated data 
first in order to gain confidence in the method. Data, reported in the literature, for 
a high rate discharge of an experimental lithium/thionyl chloride cell is used for an 
analysis. 


1 


INTRODUCTION 

This work is part of an ongoing effort to characterize the lithium/thionyl chloride 
fLi/SOC^) cell in order to better understand the physical system and improve its safety 
and performance. The Li/SOCl 2 cell is an attractive primary energy source because of 
its high energy' density (1, 2). A description of the Li/SOCU cell along with the safety 
hazards it exhibits is discussed elsewhere (2). 

A detailed mathematical model of this battery has been presented (3). This model 
is used here in conjunction with a parameter estimation technique in order to establish 
a means to estimate the electro-kinetic parameters which aid in describing the physical 
system. The assumptions and governing equations used in the model are given elsewhere 
(3). A description of the model parameters, the sensitivity of the mode] predictions 
to changes in these parameters, and estimates of these parameters based on available 
experimental data (4) is presented next. 

MODEL PARAMETERS 

The model (3) includes several reactions: the oxidation of lithium at the anode 

Li — ► Li++e“ [1] 

the reduction of SOCI 2 at the cathode 

2S0C1 2 + 4e“ — > 4C1“ + S0 2 + S [2] 

and the precipitation of LiCl on the pore surfaces of the cathode during discharge by 
the following reaction 

Li+ + Cl" ^ LiCl(s) [3] 

Due to the low solubility of LiCl in the SOCI 2 + LiAlCl* mix (5), reaction [3] is assumed 
to occur completely and instantaneously. Thus, reactions [2] and [3] can be combined to 
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yield one overall reaction representing the processes occurring within the porous cathode 

( 3 ) 

4Li + + 4e~ + 2S0C1 2 — ♦ 4LiCl + S0 2 + S [4] 


Four species are included in the model, Li + , AlCl^, SOCl 2 , and LiCl. 

The electrochemical reaction rates, depend upon the local concentrations of the 
chemical species, the potential driving force for reaction, and temperature. Analogous 
to chemical kinetics, these dependencies are not described by fundamental laws. The 
Butler- Volmer polarization equation for a flat plate electrode (6, 7) is used in the model 

( 3 ) 


* k — *o,ib,re / 




where 




[53 


Vk = $s~ $o,e - Uk.reJ 


[6] 


In Eq. [6], rjk is the potential driving force for electrochemical reaction k. The electro- 
kinetic parameters i 0 ,k,refi a a,k , and « c ,fc are the exchange current density, anodic 
transfer coefficient, and cathodic transfer coefficient for reaction k , respectively. These 
electro- kinetic parameters characterize the electrochemical reaction via Eq. [5] and are 
the parameters that are sought in this work. Note that the transfer coefficients for a 
given reaction sum to the total number of electrons involved in that reaction 


ot a ,k + <*c,* = rik [7] 

Hence, only one of the transfer coefficients need be estimated and the remaining one 
can be obtained from Eq. [7]. 
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Equation [5] must be modified slightly in order to describe the rate of reaction [4] in 
the porous carbon cathode. The rate of reaction [4] must be expressed per unit volume 
of porous electrode. This change is accomplished by premultiplying Eq. [5] by the 
available active surface area per unit volume of electrode, a. The variable a changes as 
the LiCl precipitates on the pore walls, according to reaction [4], and can be described 
by (8, 9) 



where ( is an experimentally determined parameter used to describe the morphology 
of the precipitate. In writing Eq. [8] it is assumed that all the LiCl produced from 
reaction [4] precipitates instantaneously and passivates the surface it covers. Large 
values of £ indicate needle-shaped deposits whereas small values represent flat deposits 
(9). Premultiplying Eq. [5] by Eq. [8] yields the following equation for the rate of 
reaction [4] per unit volume of porous electrode 


3 = 



(<¥)}■ '(nte)"' •>< 


RT 


V* 


) 


nfer “ p (- 2 gr» ) . 


[9] 


The effective exchange current density for reaction [4] is the product o°i 0 i re j. In all, 
nine electro-kinetic parameters. i 0 ,i,re/, a Q ,i, o Ci i, a°t 0i4)re /, o a)4 , a Ci4 . q+.i, g+ t4 , and 
q 0 ,4, and one morphology parameter, (, are used (3); two of the transfer coefficients, 
a c> i and a a< 4 , are not independent quantities because they can be calculated from a 0) i 
and a C)4 via Eq. [7]. 

Several parameters are used in the model (3) to describe the species transport. 
The diffusion coefficient, D, is used to characterize the diffusion of the salt (LiAlCl 4 ) 
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in the solvent (SOCI 2 ) due to concentration gradients in the cell. In the model (3), the 
diffusion coefficient is calculated using the following expression 

D = ft exp(3 5 /r) [10] 

where and /3 6 are constants. The transference number of charged species i, t", is used 
to describe the migration of species i in the cell due to the electric field. Two charged 
species are considered in the model, Li + and AlCl^, therefore two transference numbers 
are used, t + and ff_. Only one transference number is independent because they must 
sum to one 

*"+<■=1 [ 11 ] 

When no concentration gradients exist, the transference number of a particular species 
represents the fraction of the current carried by that species (10). When concentration 
gradients do exist, however, the transference number looses this physical interpretation 
but can be thought of as a weighting factor which determines the importance of the 
migration term in the flux expression. 

In order to use the model (3) to predict cell performance for various cell designs, 
investigate scale-up, and optimize a Li/SOCl 2 cell, estimates of the electro-kinetic 
parameters are needed. Given experimental data for a particular Li/SOCl 2 cell or half- 
cell, along with the complete cell specifications, it should be possible to use a parameter 
estimation technique to obtain these electro-kinetic parameters. A sensitivity analysis 
is conducted first to determine the relative importance of each parameter, followed 
by parameter estimation using simulated data to gain confidence in the procedure, 
and finally an experimental discharge curve is used to find the S0C1 2 electro-kinetic 
parameters. 


SENSITIVITY ANALYSIS 
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It is important to determine the sensitivity of the model predictions to changes 
in the electro-kinetic parameters. If the model predictions are relatively insensitive to 
one or more of the parameters then a fairly wide range of values for these insensitive 
parameters could be used without significantly affecting the predictions of the model. 
Such parameters need to be identified so that they can be discarded from the parameter 
estimation procedure. Also, it is often the case that the more parameters that 
are estimated, the more uncertain are the estimates due to interaction between the 
parameters, poor seeding, and round-off error (11). In addition to the identification of 
insensitive parameters, the operating conditions or situations under which the sensitivity 
to given parameters is maximum can be determined. 

The sensitivity of the model predictions to changes in parameters is determined here 
by monitoring the change in cell voltage. While holding all other parameters constant, 
the parameter of interest is perturbed slightly and the resulting change in cell voltage 
is noted. A sensitivity coefficient is calculated as follows 

_ S?-i |AE(t,)l 
m |AP ( | 

where 

A E(tj) = E(tj) - [13] 

P — 

A Pi = - V 1 - M 

■* i 

Pi and P* are the perturbed value of electro-kinetic parameter i and the reference value 
of parameter i, respectively. E(tj) is the value of the cell voltage at time tj when using 
Pi and E*(tj ) is the value of the cell voltage at time tj when using P*. In Eq. [12], n 
is the number of times over which the voltages are compared. 

Many factors come into play when calculating the sensitivity coefficients. The 
influence of a particular parameter on cell voltage varies depending on the discharge 
rate, depth of discharge, and temperature. For example, the exchange current densities, 
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*0,1, re/ and o°io,4,re/> strongly influence the initial drop in cell voltage from the open- 
circuit potential (12). On the other hand, the powers on the concentration terms, <7+.i- 
9 +) 4, and q Q )4 , have more influence later in the discharge when concentrations change 
more drastically due to localized reactant depletion. Therefore, the value of n in Eq. 
[12] is chosen so that approximately S0% of the discharge is included in the analysis. 

The results of the sensitivity analysis are shown in Fig. 1 and 2. Figure 1 shows 
the sensitivity of the cell voltage to the electro-kinetic parameters describing the lithium 
oxidation, reaction [1], over a range of discharge rates. Figure 2 is an analogous plot 
for reaction [4], The most influential parameter by far is the lithium oxidation transfer 
coefficient o a ,i. This is to be expected because a 0i i is in the exponential portion of 
Eq. [5]. The lithium oxidation exchange current density, i 0 ,\,ref > and SOCI2 reduction 
transfer coefficient are next followed by the effective SOCI2 exchange current density, 
a °*o.4,re/i and the morphology parameter, £. The electro- kinetic parameters describing 
the SOCI2 reduction display less influence than those describing the lithium oxidation. 
This could be due to the fact that the lithium oxidation occurs at the surface of a 
flat-plate electrode (a boundary condition in the model) whereas the SOCI2 reduction 
occurs throughout the porous structure of the cathode (a pseudohomogeneous region 
in the model). The overall influence of the SOCI2 reduction electro-kinetic parameters 
is, in affect, decreased because the change in voltage is distributed over the thickness of 
the cathode as opposed to the change occurring all at one point. The exponents on the 
power terms, 9+4, 9+4, and q 0yi show little influence at the discharge rates investigated. 
They are discarded from most of the parameter estimation work presented next; they 
are set to the corresponding stoichiometric value as is often done. 

The sensitivity of model predictions to changes in two of the transport parameters, 
and /?6, was briefly investigated as shown in Fig. 3. At high current densities, the 
transference number shows substantial influence and the effects of diffusion become 
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important. However, the sensitivity to these parameters is reduced as the current is 
reduced because the electrochemical reaction rates become controlling. This result 
suggests that a parameter estimation procedure aimed at finding the electro-kinetic 
parameters should be performed with low-rate data so that species transport will not 
interfere with the accuracy of the estimates. 

PARAMETER ESTIMATION 

Klinedinst and Domeniconi (4) investigated the performance of a single Li/SOCl 2 cell 
with various carbon electrode thicknesses and operated over various discharge rates. 
Their experimental test cell has some features which make it attractive for obtaining 
experimental cell voltages for comparison with the model predictions here. The sin- 
gle cell, corresponding to the single cell model (3), is immersed in excess electrolyte 
and hence the cell temperature remains nearly constant throughout the discharge. The 
electrolyte is introduced into the test cell shortly before the beginning of the test (ap- 
proximately 10 minutes before) suggesting that little or no secondary LiCl film is present 
on the Li electrode. Also, the anode, separator, and cathode are pressed together with 
a weight so that no reservoir exists between the separator and cathode. 

The cell model (3) was changed to simulate this test cell for the parameter 
estimation. The film and reservoir regions were eliminated from the model by setting 
the film porosity equal to the separator porosity, removing the reservoir calculations 
from the model, and modifying slightly the boundary conditions at the reservoir /porous 
electrode interface to be those for the separator /porous electrode interface. The model 
inputs for simulating the test cell of Klinedinst and Domeniconi are listed in Table I. 

Initially, simulated data for both a 1 mA/cm 2 and a 31 mA/cm 2 discharge were 
used for the parameter estimation in order to gain confidence in the procedure and guide 
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future experimentation. Experimental data, for a 31 mA/cm 2 discharge, were used to 
estimate the SOCl? reduction electro-kinetic parameters for this high-rate discharge. 
Reformulation of the Butler- Volmer Equation: 

To reduce the interaction between parameters in the least -squares procedure, it 
is sometimes necessary to reformulate the model equations. Parameter interaction is 
evident when small changes in the parameters produce similar changes in the objective 
function. Interaction often occurs when the parameters appear as a product somewhere 
in the governing equation. Interaction can be reduced if the equation is rewritten so 
that the parameters appear in separate terms. 

As shown in Eq. [9], the electro-kinetic parameters in the Butler- Volmer equation 
appear in product terms as opposed to separate terms. The Butler- Volmer equation 
can be rearranged by first taking the natural logarithm and then taking the exponent 
of the various terms (12). For the parameter estimation used here, Eq. [9] is written as 
follows 

j = exp [/i + h] - exp [/i + / 3 + /4 + /s] - 


exp [/i + /2 + /e] + exp [/j + fs + h + h + /e] [ 15 ] 


where 


/i = ln(a°i 0)4 , re/ ) [ 16 ] 

h = - *2 - U,,r'f) [17] 

h = ?+, 4 In(-T-) [18] 

C re/ 

f 4 = q a , 4 ln(— ^ — ) [19] 

c o, ref 

h = ( *' - ** ~ Pol 

/• = £ ln ( 5 V i ) l 21 l 
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Equation [5] for reaction [1] is rewritten in a similar fashion. As a result of this refor- 
mulation, the parameters used in the parameter estimation procedure are ln(a°i 0 , 4 ire ^). 
^c .4 1 ^a,i 5 and ln(i 0( i ire /). 

Introduction of Noise into Simulated Data: 

In order to test the ability of the parameter estimation procedure to predict the 
model parameters within reasonable accuracy, random noise was added to the simulated 
data. A random number generator from the IMSL library (13 < was used to generate a 
normally distributed number, <p, with mean zero and variance one (i.e., ~N(0,1)). The 
noise was added to the predicted value of the cell potential at time tj, E(tj), as follows 

E{tj ) with noise = E(tj^) + ipZ [22] 

where Z is the average level of noise (in mV) presumably inherent in the measurement 
of the cell voltage. Three levels of Z were investigated here, 0.5 mV, 1 mV, and 2 mV. 
It is realistic to assume that the experimental measurements of the cell voltage would 
be accurate to within 2.0 mV. 

Calculation of Confidence Intervals for Parameter Estimates: 

An estimate of a parameter has little meaning unless it is accompanied by some 
approximation of the possible error it possesses. Here, corndence intervals for the 
parameter estimates can be estimated by using the following formula (14) 

Pi = Pi ± sp. [23] 

In Eq. [23], Pi is the estimate of P,-, <( 1 _i) > ( n _ m ) is the t-distribution at the (1 — ^) x 100 
percent confidence level, (n — m) is the degrees of freedom, n is the number of 
observations, m is the number of parameters estimated, and #p t is the estimate of the 
variance of P, which is calculated from the sum of squares of the error, as follows: 

s p i — \J C u s E 


[24] 
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In Eq. [24], Cu is an element of the matrix (J-Tj ) -1 (i.e., inverse of the approximate 
Hessian, see appendix B in ref. 15). The estimate of the variance in the error, s 2 E . is 
calculated from the sum of squares of the residual errors ( SSR ) as follows 

2 SSR Ylj=\(E(tj)obs ~ E(tj)pred ) 2 , , 

*E = = — L 2 ^> 

n — m n — m 

Parameter Estimation Using High-Rate Simulated Data: 

The use of high-rate (e.g., 30 mA/cm 2 ) discharge data for parameter estimation 
seems to offer several advantages. High rate discharges occur over a relatively small 
time as opposed to low rate discharges (e.g., about 1 hour for a 30 mA/cm 2 discharge as 
opposed to approximately 50 hours for a 1 mA/cm 2 discharge of a Li/SOCl 2 cell). This 
means less time in the lab and, perhaps, less computer time to simulate the discharge. 
Often, the shapes of voltage-time curves for high-rate discharges are much more curved 
than the relatively flat voltage-time curves for low-rate discharges thus offering more 
information for the parameter estimation. 

Here, however, the use of high-rate data exhibited several disadvantages. First, 
the influence of species transport becomes important at the higher discharge rates. For 
example, at 100 mA/cm 2 the sensitivity coefficient, as calculated using Eq. [12], for the 
transference number is 0.3 which is second only to the sensitivity coefficient of a a 4 , 0.53. 
Due to the uncertainty in the transference number and the lack of data for the diffusion 
coefficient (see appendix A in ref. 15) it is best to stay clear of those operating regions 
where species transport is a controlling aspect of the cell discharge. Secondly, at high 
rates much interaction exists between the electro-kinetic parameters. This interaction 
can be seen in Fig. 4 which is a plot of the time-dependent sensitivity coefficients 
(i.e., the elements of the Jacobian matrix, as shown in appendix B in ref. 15) for a 
31 mA/cm 2 simulated discharge. These sensitivity coefficients are calculated using Eq. 
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[12] except that the perturbation in each parameter, A P, is not scaled by P* as in Eq. 
[14] but is calculated as follows 

A P.; = Pi - P * [26] 

Note that three of the five parameters, ln(a 0 ? 0 , 4 ,re/)i a a .j, and ln(i 0] i )r£ /), influence the 
cell voltage in the same fashion. This similar influence (i.e., interaction) is shown by the 
three straight parallel lines (14) in the figure. Thus, only one of these three parameters 
can be determined independently. Note also that the two lithium oxidation parameters 
cannot be individually obtained at high rates. 

Another indication of this interaction is the condition number of the approximate 
Hessian (11) (see appendix B in ref. 15). The condition number of a matrix is the ratio 
of its largest eigenvalue to its smallest eigenvalue. A matrix is called ill-conditioned when 
its condition number is large. The larger the condition number, the more substantial 
the change in the solution even with very small changes in the data. In essence, a 
large condition number here means that a unique set of electro-kinetic parameters is 
more difficult to find; many sets will reduce the objective function to nearly the same 
value. The condition number for the approximate Hessian here was on the order of 
1 x 10 18 . The reason for this ill-conditioning could be the fact that transport becomes 
more important at higher rates which serves to decrease the influence of changes in the 
electro-kinetic parameters on cell performance. 

Though interaction between some of the parameters exists at high rates, several 
attractive features can be noted. If the lithium oxidation parameters, i 0 ,i,ref and o Qj i , 
are known a priori via independent means, then the SOCI 2 electro-kinetic parameters 
can be determined using high-rate data. This ability is shown in Table II where 
ln(a°t 0) 4 )re /), a c> 4 , and f are estimated from simulated 31 mA/cm 2 data with 0.5 
mV added random noise. This ability is important because the SOCI 2 electro-kinetic 
parameters will probably change with discharge rate due to changes in the active surface 
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area in the porous electrode as discharge proceeds. In addition, the sensitivity of the 
predicted cell voltage to changes in £ seems to be greater at high rates than at low rates 
thus enabling better estimates of f. This can be seen by comparing the sizes of the 
confidence intervals for £ here with those in the next section. 

Parameter Estimation Using Low-Rate Simulated Data: 

Parameter estimation on low-rate data (e.g., 1 mA/cm 2 ) remedies some of the 
problems encountered with the high-rate data. The interaction between parameters 
is decreased. Figure 5 is a plot of the sensitivity coefficients versus time for the 
1 mA/cm 2 simulated discharge. The curves have differing slopes and shapes, as 
opposed to those for the high-rate case as shown in Fig. 4, indicating that the various 
parameters affect cell voltage differently and hence can be determined independently. 
Note that these coefficients are calculated using the known values of the electro-kinetic 
parameters. Table III lists the condition numbers for the approximate Hessians for 
various combinations of the parameters; the various combinations are indicated by case. 
The condition numbers for cases A and E are low suggesting that parameter estimations 
for either of these cases ought to be straightforward and the confidence intervals on the 
parameter estimates small (11). Here, the term straightforward parameter estimation 
means that the search proceeds to the solution given a variety of initial guesses that 
can be “far” from the correct solution. More difficult parameter estimation runs involve 
multiple initial guesses in order to make sure the estimates obtained with a particular 
initial guess are not those at a local minimum. A comparison of cases A and B shows that 
estimates for all three SOCI 2 electro-kinetic parameters will be more difficult to obtain 
than estimates for ln(a 0 i 0) 4 ire /) and a Ci4 alone. Case G shows that the eight-parameter 
problem is extremely ill-conditioned, meaning much interaction between parameters and 
very large confidence intervals on the parameter estimates. Therefore, q+,i, q+, 4. and 
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q 0 , 4 are omitted from the following parameter estimation work and are set to their 
stoichiometric values as in Table II of ref. 3. 

The simulated data for the 1 mA/cm 2 discharge investigated here, with added noise, 
are shown in Fig. 6 . The model inputs are listed in Table I. The data in Fig. 6 spans 
approximately 50% of the discharge. Although predicted discharge times are much 
longer for low-rate than for high-rate, the solution times for low-rate cases may actually 
be smedler than those for high-rate cases. Solution times depend upon the number of 
spatial nodes, the size of the time steps, and the number of iterations (at each time) to 
converge upon the solution. At low-rates, changes in the dependent variables (electrolyte 
concentration, current, and potential) in space and time occur much less rapidly than 
at high-rate. Therefore, fewer node points and larger time steps can be used at low- 
rate than at high-rate to maintain three-digit accuracy. For example, 81 nodes and a 
maximum time step of 37.5 seconds were used to maintain three-digit accuracy for the 
31 mA/cm 2 discharge results reported earlier. Here, for the 1 mA/cm 2 discharge, 21 
nodes and a maximum time step of 2000 seconds were sufficient. This result eliminates 
the advantage of using high-rate data to reduce solution time; the only disadvantage to 
using low-rate data is the length of time it takes to collect the experimental data. 

The results of the parameter estimation runs are shown in Tables IV - VII. As 
expected, the runs for cases A and E were straightforward with narrow confidence 
intervals for the parameter estimates. As shown in Tables IV and V, the initial guesses 
are at different extremes and up to two orders of magnitude away from the true solution. 
Notice how the confidence intervals increase as the noise level increases which is to be 
expected. Case E may be useful when estimates of the SOCI 2 electro-kinetic parameters 
have been obtained independently. Similarly, Case A may be needed when estimates 
of the lithium oxidation electro-kinetic parameters are available and an estimate of £ is 
known from independent studies. Table VI lists parameter estimation runs for case B, all 
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three SOCI 2 reduction electro-kinetic parameters. This case, as shown in the previous 
section, can be used to determine both low-rate and high-rate values of a°i o i re f , a Ci 4 , 
and f when estimates of a 0t i and ln(i 0) i (re j) are known. Notice the increase in the size 
of the confidence intervals from Table V to Table VI. The estimation of the additional 
parameter, £, adds greater uncertainty to the parameter estimates, as expected. Also, 
judging from the large increase in condition number from case A to case B in Table III, 
some interaction between these parameters is evident. 

Table VII lists parameter estimation runs for case D. The ability to estimate all five 
electro-kinetic parameters at once, as in case D, is attractive. Simultaneous estimation 
of all parameters in the model usually provides parameter estimates possessing superior 
predictive capability. Often, independently determined parameter estimates do not 
predict the behavior of the entire unit very well. For example, independent estimates of 
the lithium oxidation parameters and SOCI 2 reduction parameters may not be useful 
for predicting the behavior of the Li/S0Cl2 cell. These independent estimates might be 
obtained from rotating disk electrode experiments and a corresponding mathematical 
model. In Table VII, note that for the data with added 0.5 mV noise the confidence 
intervals on certain parameter estimates are unacceptably large. This uncertainty in the 
predictions is due to parameter interaction (note the relatively large condition number 
in Table III) even at this low rate of discharge. The five parameters can be obtained 
with no noise, as shown in Table VII, but it is unrealistic to expect exact measurements 
and the model to be an exact description of the physical system. It may be useful to 
obtain all five parameters simultaneously, even though the confidence intervals are large, 
because the estimates should still provide an excellent fit of the model predictions to the 
experimental data (notice the small SSR in Table VII). The parameter estimation runs 
for cases B and D were not straightforward. Certain initial guesses of the parameters 
detoured the search algorithm to local minimums. For example, for the case D runs, 
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the initial guesses of the exchange current densities needed to be within an order of 
magnitude and the transfer coefficient guesses needed to be within about half of their 
true value in order to reach convergence. The algorithm could not escape these local 
minimums within a reasonable number of iterations (e.g., 60 iterations). 

An alternative to predicting all five parameters simultaneously would be to obtain 
estimates for some of the parameters independently (for example, the lithium oxidation 
parameters) and then use these estimates in the cell model to obtain the remaining 
parameters (the SOCI 2 reduction parameters) using data from the entire cell. In 
addition, separate data for each electrode in the Li/SOCU cell can be obtained by 
constructing half-cells. For example, a lithium electrode half-cell could be constructed 
from a Li/SOCl 2 cell by placing a reference electrode near the surface of the lithium 
electrode. The lithium electrode is the working electrode, a second lithium electrode 
(with high impedance) placed near the surface of the working electrode could act as 
the reference electrode, and the porous carbon cathode is the counter electrode. The 
potential of the lithium electrode, as it changes with time, could then be measured 
relative to the reference electrode. Predictions of the Li electrode potential could be 
obtained from a model of the working electrode / reference electrode system, consisting 
of equations presented earlier for the entire Li/SOCb cell. These predictions could 
then be compared with the experimental values to obtain estimates of the electro-kinetic 
parameters of the lithium electrode. 

Parameter Estimation Using Available Experimental Data: 

The analyses using simulated data in the previous sections have shown the need 
to collect low-rate data from a well defined experimental cell. Unfortunately, no such 
experimental data is currently available. Various experimental voltage-time curves have 
been reported (16, 17, 18, 19) but it is difficult to use this data due to the lack of 
cell specifications given (e.g., porous electrode porosity and thickness, initial electrolyte 
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concentration, temperatures, and separator porosity and thickness). In addition, large 
differences exist in the experimental data which can be attributed to the different 
electrode materials, initial electrolyte concentrations, operating conditions, and cell 
dimensions used. 

Using data from Klinedins* (20), for a 31 raA/cm 2 discharge of an experimental cell 
(4), estimates of the SOCU reduction electro-kinetic parameters were obtained. The 
lithium oxidation electro-kinetic parameters were arbitrarily set to the values listed in 
Table I of ref. 3. These values are well within ranges commonly used in the literature 
for electrochemical reactions. Initially, a by-hand trial and error search was conducted 
to find electro-kinetic parameters giving model predictions close to the experimental 
results. It was hoped that this search would provide initial guesses in the vicinity of 
their actual values to aid in the parameter estimation. Two sets of initial guesses were 
tried, as shown in Table VIII as Run # 1 and Run # 2, and the same final parameter 
values were obtained to within the confidence limits shown. The values are listed in 
Table VIII and the predicted voltage curve, obtained using the estimates from Run # 1, 
is superimposed on the experimental data in Fig. 7. Note that the confidence intervals 
for a 0 i O) 4 ) re/ and ( in Run # 1 and Run # 2 are large relative to the estimated values 
of the parameters. To reduce the uncertainty in the estimate of a°t 0 , 4 ,re/, s was set to 
a value of 0.02 (a value between the estimates of £ in Run # 1 and Run # 2 ) for Run # 
3. This procedure (estimating all three SOCI 2 electro-kinetic parameters, setting £ to 
an average value, and then estimating a°t 0 ,4,re/ and a c ,4) gives estimates of a°i 0 ^ tre f 
and a C) 4 with reasonable confidence intervals and an estimate of £ in the vicinity of 
its best value. 

Notice in Fig. 7 that the first experimental data point used occurs some time 
after the beginning of the discharge. This is due to the fact that the beginning behavior 
of experimental discharge curves are very different. Some exhibit a very sudden and 
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large initial decrease followed by a time of voltage recovery and is known as the voltage- 
delay phenomena. Others show a steady decrease, with an initial large decrease, as 
the model predictions do. The voltage-delay phenomena in Li/S0Cl2 cells has been 
studied by Dey (21) and has been associated with the LiCl film. Slow voltage recovery 
is associated with slow mechanical disruption of the film. The model does not account 
for this voltage-delay phenomena and therefore initial experimental data were excluded 
from the parameter estimation. 


CONCLUSIONS 

The sensitivity analysis and parameter estimation performed here for the model 
presented by Evans et al. (3) yields several conclusions. The test cell of Klinedinst and 
Domeniconi (4) offers a constant temperature environment, a single cell (as opposed 
to a stack), elimination of film and reservoir regions, and well-defined thicknesses of 
the cell components. Low-rate data should be used to obtain estimates of the electro- 
kinetic parameters so that interaction between parameters is small. Low rate data 
is just as efficient to use in the parameter estimation procedure as is high rate date 
because the number of nodes can be decreased and the time step increased so that the 
solution times are comparable. Also, at low rates the effects of species transport are 
minimal. Once estimates of the lithium oxidation parameters are available, the electro- 
kinetic parameters describing the SOCI 2 reduction on the porous carbon electrode 
can be determined even for high discharge rates. These parameters will change with 
electrolyte composition, type of carbon used for the electrode, and temperature. 
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NOMENCLATURE 

a specific electroactive surface area of the porous electrode, m -1 

a 0 initial value of a. m _1 

c, concentration of species i, mol/m 3 

Cj ntt initial concentration of electrolyte, mol/m 3 
c i, re f reference concentration of species i, mol/m 3 
Ca diagonal element of approximate Hessian 

D diffusion coefficient of the binary electrolyte, m 2 /s 

E cell voltage. V 

F Faraday’s constant, 96487 C/mol of electrons 

I* current density due to electrochemical reaction it. A/m 2 

iok,rtf exchange current density of electrochemical reaction k at c rt /, A/m 2 

j current transferred between phases, A/m 3 

n number of data points 

n* number of electrons transferred in the electrochemical reaction k 

m number of parameters 

Pi parameter i 

qik cathodic reaction order of species i in reaction k 

R universal gas law constant, 8.314 J/mol-K 

se estimate of the variance of the residual error 

sp. estimate of the variance in the estimate of Pj 

Si thickness of region 1, m 

SCi sensitivity coefficient for parameter i 

SSR sum of squares of the residual error 

T cell temperature, K 

tj time j , sec 

Af mn x maximum time step, sec 

<7 transference number of species i 

Uk,ref open circuit potential of reaction k based on the reference 
concentrations, V 

Greek 

ot ak anodic transfer coefficient for reaction k 

a c k cathodic transfer coefficient for reaction k 

e porosity of porous material 

t° initial porosity of the porous electrode 

T)k overpotential for reaction k at electrode surface (V — — C'k,ref), V 

if morphology parameter 

$ solution potential, V 

$o solution potential at electrode surface, V 
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Table I. Model inputs used to simulate a 1 mA/cm 2 discharge of an experimental 
test cell (4). 


Numerical parameters: 

Total number of node points = 21 

Maximum time step (A t maz ) — 2000 seconds 

Physical properties and initial conditions: 

Cell temperature ( T ) = 25°C 

Initial electrolyte concentration (c tntt ) = l.S M 

Reference electrolyte concentration ( c re j ) = l.S M 

Transference number (t* ) = 0.9 


All other values as in Tables I and II of ref. 3. 
Cell dimensions: 


Thickness of cathode (5 pe ) 

rr 

0.05 cm 

Initial porosity of cathode (e°) 

=r 

0.S45 

Thickness of separator ( S s ) 

=r 

0.0127 cm 

Porosity of separator (e a ) 

= 

0.7 

Thickness of reservoir (S r ) 

= 

0 cm 

Thickness of film ( S/) 

= 

0 cm 

Lithium oxidation electro-kinetic parameters: 

i o,l, ref 


0.002 mA/cm 

<*a, 1 

= 

0.5 

a c,l ( = 1 - Ofl . l ) 

= 

0.5 

Q+,1 

= 

1.0 

SOCI 2 reduction electro-kinetic parameters: 

a °*o,4,re/ 

= 

0.005 mA/cm : 

<*c,4 

= 

1.0 

°a, 4 ( ~ 2 - 0 Ct 4 ) 

= 

1.0 

t 

= 

0.05 

q+A 

= 

1.0 

Qo, 4 

ss 

0.5 


Table II. Estimated values of the SOCI2 reduction electro- kinetic parameters 
using high rate (31 mA/cm 2 ) simulated data with added noise. 


Actual Values*: 

o°?o,4 i re/ = 0.0441 A/cm 3 
a C) 4 = 0.6846 
£= 0.05 


Z (mV) 

Starting Values 

Estimated Values 
(95% confidence, 7 = 0.05) 

0.5 

ln(a°? 0i 4,re/)= -5.423 

-3.073 ± .0609 


(a°i 0 ,4,ref = 0.00441) 

( C\ (Ufii +0.00291 \ 
^O.U463 —0.00274 J 


Ot c ,4 = 0-171 

0.681 ± 0.0361 


0.01 

0.0487 ± 0.01356 


The lithium oxidation electro-kinetic parameters were held fixed here at 
the following values: 

*0,1 , ref = 1.58 x 10 -3 



Table III. Condition numbers for the approximate Hessians of the least squares 
objective function for the 1 mA/cm 2 simulated discharge. Large condition numbers 


indicate parameter interaction. 





Case 

Approximate Hessian of 

Condition 

Number 


A 

ln(a°i 3 a, ref) and a c , 4 

5.62 x 10 2 

B 

ln(a 0 i O) 4 ire /), a Ci4 , and £ 

1.18 x 10 5 

C 

ln(°° l M,re/), <* C) 4, $, and «a,l 

1.14 x 10 6 

D 

ln(a l o, 4 ,re/)i a c,Ai ^a,ii and ln(i 0 i )re ^) 

6.16 x 10 7 

E 

®a, i and ln(z 0i i >re /) 

4.30 x 10 1 

F 

ln(fl°?c),4,re/)? a c,4i a a,li and ln(i 0) i )rc /) 

1.03 x 10 6 

G 

ln(Q°lo,4,re/)i a c, 4) O t a,li ln(io,l,re/)5 9+,l? 

<Z+, 4, and g 0) 4 

2.45 x 10 12 



Table IV. Estimated values of the lithium oxidation electro-kinetic parameters 


using simulated 1 mA/cm 2 data with added noise. 



Actual Values: 

io,i,r e / = 0.002 A/cm 2 
<* a ,i = 0.5 

(Case E) 

Z (mV) 

Starting Values 

Estimated Values 
(95% confidence, 7 = 0.05) 

0.5 

ln(i 0 > i t re/)= -S.517 

-6.223 ± .062 


(*o,l t re/ = 2.0 XlO 4 ) 

( 0.00198 


1= 0.15 

0.508 ± 0.072 

1.0 

ln(t 0ili re/)= -S.517 

-6.144 ±0.105 


(*o,l, re / = 2.0 XlO -4 ) 

fn nn 9 i 4 ^ +2.375 x 10 -4 \ 
^0.002145 _2.139 X io -4 J 


a a,l= 0.15 

0.432 ± 0.126 

1.0 

lri(*o,l,re/) = -10.82 

-6.144 ± 0.100 


(*o,l,re / = 2.0 XlO -5 ) 

(0 002146 


o a>1 = 0.85 

0.432 ±0.118 

2.0 

ln(l 0 ,i,re/)= -8.517 

-6.263 ± 0.211 


(*o,l,re/ = 2.0 XlO 4 ) 

(o.OO 1906 ^ 3^625 x 10 “* ) 


a a ,i= 0.15 

0.5337 ± 0.248 




Table V. Estimated values of the exchange current density and transfer coefficient 
of the SOCI 2 reduction reaction using simulated 1 mA/cm 2 data with added noise. 


Actual Values: 

a °i 0 ,4,ref = 0.005 A/cm 3 

d' Ci 4 = 1.0 

(Case A) 

Z (mV) 

Starting Values 

Estimated Values 
(95% confidence, 7 = 0.05) 

0.5 

ln(a°i 0 i 4 >re /)= -2.996 

-5.293 ± 0.0732 


(a°io,4,ref = 0.05) 

(o.00503 


ft Ci 4 = 0.25 

0.9981 ± 0.0212 

0.5 

ln(a°i 0 i 4 ire /)= -9.903 

-5.293 ± 0.0732 


( a°io, 4 ,ref = 5.0 xlO -5 ) 

( 0.00503 ii.11 ^ f o z : ) 


ct c , 4 — 1-75 

0.9981 ± 0.0212 

1.0 

ln(a°i 0 i 4 ire /)= -0.6931 

-5.239 ± 0.123 


(0 io,4,ref — 0.5) 

( n onto +6.948 xlO " 4 \ 
^0.0053 _6 i 44 x 10 -4 J 


a C) 4 = 0.25 

0.9842 ± 0.0357 

2.0 

ln(a°z 0 ) 4 )re y)= -9.903 

-5.212 ± 0.234 


(a°io,4,reJ = 5.0 XlO" 5 ) 

( 0.00545 ± 8 : 88114 ) 


a C)4 = 1.75 

0.9727 ± 0.0677 




Table VI. Estimated values of the SOCI 2 reduction electro-kinetic parameters 
using simulated 1 mA/cm 2 data with added noise. 


Actual Values: 

Q-° i 0 ,4 ,re f = 0.005 A/cm 3 
a c 4 = 1.0 
£= 0.05 

(Case B) 

Z (mV) 

Starting Values 

Estimated Values 
(95% confidence, 7 = 0.05) 

0.5 

ln(a°io, 4 ,re/)= -7.601 

-5.289 ± .129 


(a°io,4,ref = 5.0 xlO -4 ) 

( n +6.95 xlO” 4 \ 

^0.00505 _ 611xl0 -<j 


o Ci 4 = 0.25 

1.001 ± 0.075 


£= 0.01 

0.0494 ± 0.0184 

1.0 

ln(<z 0 i 0 i 4 )re /)= -9.903 

-5.257 ± .189 


(a°i 0 ,4,ref = 5.0 xlO -5 ) 

( n nn^oi ++085 X 10- 3 \ 
^0.00521 _ g 9g x 10 -4 J 


p 

0 

II 

01 

0.9707 ± 0.1238 


£= 0.01 

0.0535 ± 0.0313 



Table VII. Estimated values of all five electro-kinetic parameters describing the 
Li/S0Cl2 cell using simulated 1 mA/cm 2 data with added noise. 



Actual Values: 

a°i 0 A,ref = 0.005 A/cm 3 
a c , 4 = 1.0 
£= 0.05 

*0,1 , ref = 0.002 A/cm 2 
c* a ,i = 0.5 

(Case D) 

Z (mV) 

Starting Values 

Estimated Values 
( 95 % confidence, 7 = 0.05) 

0.0 

ln(a°i 0)4 , re/ )= -7.601 

-5.298 


( a *o,4 , ref — 5.0 X 10 4 ) 

(5.0 x 10" 4 ) 


« c ,4= 0.5 

1.00 


(= 0.025 

0.05 


kl(lo,l,re/) :=: -&.517 

-6.215 


0o,l,re/ = 2.0 X 10 4 ) 

(2.0 x 10~ 4 ) 


a a ,i= 0.25 

0.5 


SSR = 1.0 x 10~ 18 


0.5 

ln(a 0 i o ,4,r e /)= -7.601 

-5.031 ± 2.078 


( a°io,4,ref = O.o X 10“ 4 ) 

^0.00653 Jomn) 


a C 4= 0.5 

1.00 ± 0.253 


0.025 

0.06-56 ± 0.0678 


^( l o,l,re/) = -S.517 

-6.523 ± 1.292 


(*o,l,re/ = 2.0 X 10 4 ) 

( n "t~0. 00388^ 

^0.00147 —0.00107 J 


C*a,l = 0.25 

0.374 ± 1.80 


SSR = 1.067 x 10~ 5 



Table VIII. Estimated values of the SOCI 2 reduction electro-kinetic parameters 
of the Li/SOCl 2 cell using available experimental data at 31 mA/cm 2 . 


Starting Values 


Estimated Values 
(95% confidence, = 0.05) 


Run # 1 (Case B) 


0.95 x 10~ 4 


ln(a°t 0) 4 , r< /) — -3.0 
(fl 0 t'o,4.re / = 0.05) 

a c , 4 = 0.5 
t= 0.05 


-1.51 ±5.78 


0.221 


+71.3 


0.657 ±0.134 
0.0112 ±0.0724 


Run # 2 (Case B) 


0.96 x 10- 


ln(a 0 t O) 4 . r< /)= -7.601 
(a°i 0 ,4.ref = 5.0 x 10" 4 ) 
a C)4 = 1.5 

£= 0.1 


-2.3 ±1.95 

ion ±0.604 
I 0.100 _q.086 

0.64 ± 0.139 
0.028 ± 0.0736 


Run # 3 (Case E) 


0.95 x 10“ 4 


ln(a 0 i Oi 4 , re/ )= -3.0 


( a *o,4 


0.05) 


-2.00 ± 0.655 

( o ±0.125 

(0.135 _q 065 

0.646 ± 0.094 


a Cl 4= 0.5 
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electro-kinetic parameters characterizing the lithium oxidation, 
reaction [1]. 

Sensitivity of the predicted cell voltage to changes in the 
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discharge evaluated at the exact values of the electro-kinetic parameters. 
Parallel lines indicate parameter interaction. 
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A simulated 1 mA/cm 2 discharge of a Li/SOCL cell with three levels 
of added random noise. 

A comparison of model predictions to experimental data at 31 mA/cm 2 . 
Electro-kinetic parameters are as shown in Table VIII, Run # 1. 



0.6 


- 0.1 


t /v-3 ^ &, , -JL VJ Ify 



30.0 40.0 50.0 60.0 

i n (mA/cm 2 ) 
















APPENDIX B 



A Thermal Analysis of a Spirally Wound Battery Using a 
Simple Mathematical Model 


by 


T. I. Evans* and R. E. White** 
Department of Chemical Engineering 
Texas A&M University 
College Station, Texas 77843 


submitted as a technical paper 
to the 
Editor 

Journal of the Electrochemical Society 
10 South Main Street 
Pennington, New Jersey 08534 

April, 1988 


Key Words: discharge, energy balance, lithium/thionyl chloride, spiral, thermal 
runaway 

* Electrochemical Society Student Member 
^^Electrochemical Society Active Member 



ABSTRACT 


ORIGINAL PAGE IS 
OF POOR QUALITY 


A two-dimensional thermal model for spirally wound batteries has been de- 
veloped. The governing equation of the model is the energy balance. Convective 
and insulated boundary conditions are used and the equations are solved using a 
finite element code called TOPAZ2D. The finite element mesh is generated using 
a preprocessor to TOPAZ2D called MAZE. The model is used to estimate tem- 
perature profiles within a spirally wound D-size cell. The model is applied to the 
lithium/thionyl chloride cell because of the thermal management problems that this 
cell exhibits. Simplified one-dimensional models are presented which can be used to 
predict best and worst temperature profiles. The two-dimensional model is used to 
predict the regions of maximum temperature within the spirally wound cell. Normal 
discharge as well as thermal runaway conditions are investigated. 



INTRODUCTION 


Many battery systems produce heat as they are discharged due to the exother- 
mic nature of the electrochemical reactions occurring. This build up of heat can 
result in dangerous conditions depending upon the cell design and physical system. 
Often, thermal modeling of these battery system; is accomplished by assuming that 
the different components of the battery interior can be represented by a homoge- 
neous region (i.e. core region) having average properties. This approach is useful 
because it greatly simplifies the analysis, however, it may not provide adequate 
results when the different components of the battery have vastly different thermal 
properties. Often, the material used for separating the anode from the cathode has 
a much lower thermal conductivity than the electrodes themselves. This condition 
may lead to heat management problems in spirally wound designs when the heat 
flow is primarily in the radial direction. Heat management could be improved in 
the spiral design if heat flow can proceed out the spiral along the more conductive 
materials; however, the spiral path offers a much greater distance to traverse as 
opposed to the more direct radial path. This work shows that conduction of heat 
out the spiral, during normal operation, reduces cell temperatures but does not 
substantially reduce the temperature drop in the cell. However, when hot spots are 
present in the cell conduction out the spiral is substantial and the temperature drop 
in the cell is reduced significantly. 

The Li/SOCl 2 cell is an example of a spirally wound D-size cell that has a 
thermal management problem. It is an attractive primary energy source because 
of its high energy density (1, 2). However, researchers have observed that high 
discharge rates and high temperatures promote thermal runaway in these cells (1- 
5). Venting of toxic gases and explosions, due to the rapid pressure build up at 
higher temperatures, have been reported (1-5). 
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The purpose of this paper is to present a model that can be used to predict 
temperature profiles, in time and two spatial dimensions, within a spirally wound 
D-size Li/SOCl 2 cell and to estimate the ability of this design to conduct heat out 
of the cell. One-dimensional approximations to the two-dimensional analysis will be 
investigated for comparison purposes. Predictions from a one-dimensional model 
for a C-size cell are compared with experimental data (6) to learr. more about 
the thermal behavior of Li/S0Cl2 cells. Finally, hypothetical thermal runaway 
situations are simulated using the two-dimensioned model and a one-dimensional 
approximation in order to determine if the spiral design offers any thermal advantage 
under such conditions. 

Description of a spirally wound electrochemical cell 

Figures 1 and 2 show two different views of a spirally wound electrochemical 
cell. The cell is constructed by inserting the cell roll into a cylindrical can as shown 
in Fig. 1. Electrolyte is poured into the can and fills the porous regions of the 
roll so that a reservoir of electrolyte remains at the top of the roll. Electrodes are 
constructed by impregnating the electrode materials on the both sides of the current 
collector mesh (usually a nickel mesh in the Li/SOCl 2 cell). The anode current 
collector is usually spot welded to the side of the can and tabs connected to the 
cathode current collector are welded to the center post protruding through the can 
lid. The cell roll is constructed by stacking the anode, separator, and cathode on 
top of each other and then rolling up the stack. A top view of the roll is shown 
in Fig. 2 - which is the two-dimensional plane of interest. The current collectors, 
embedded in the middle of each electrode, are not shown in Fig. 1 or 2. 

Description of a Li/SOCk cell 
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In a typical Li/S0Cl2 cell, the anode is Lithium, the cathode is a porous carbon 
structure, and the separator is a glass matting material. The electrolyte consists 
of lithium tetrachloroaluminate (LiAlCl 4 ) in thionyl chloride (SOCI 2 ). During 
discharge, Li is oxidized 

Li — ♦ Li + -f e~ [T 

and SOCI2 is reduced 

2S0C1 2 + 4e" — ♦ 4C1“ + S0 2 + S [2’ 

The SOCI 2 is the active material at the cathode; the porous carbon structure 
provides the sites upon which the SOCI 2 is reduced. This redox couple generates 
heat as the cell is discharged and this heat must be conducted away at a fast enough 
rate so that unsafe conditions will not develop. 

This brief description is a simple analysis of the actual physical system. It has 
been observed that the Li/S0Cl2 cell is a complex chemical and electrochemical 
system involving an unknown number of reactions (2). Szpak and Venkatsetty (7) 
state that modeling the Li/S 0 Cl 2 cell is difficult due to the complexity of this 
physical system. They observe that as the cell discharges, the temperature and 
pressure of the cell change, the volume of electrolyte decreases, the electrolyte 
composition varies, and new phases appear. However, much can be learned by 
modeling a simplified picture of the actual system. 

Previous models 

Previous workers (8, 9) have presented thermal models for battery systems. 
Bernardi et al. (8) give an excellent development of a general energy balance for 
electrochemical systems. They present a complete energy balance which requires the 
knowledge of concentration profiles and current fractions within the cell. However, 
the energy balance is formulated by assuming that the temperature of the battery is 
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uniform and changes only with time. That is, flat temperature profiles are assumed 
a priori. Lee ei al. (9) present a three-dimensional model for battery systems 
involving a number of batteries in one unit. However, they treat the interior 
region of each battery as a homogeneous phase having effective average properties. 
Therefore the model cannot be used to predict the effects of the arrangement of the 
cell components on the temperature distribution. 

Several models (10-19) have been presented which specifically address the 
thermal behavior of Li/SOCl 2 cells. Parnell and Szpak (10) present a thermal 
model for thin, disc-shaped Li/SOCb cells in one dimension. They consider heat 
generation due to the polarization of the electrodes and interna] cell resistance. 
A differential energy balance, which includes various terms representing different 
contributions to the heat generation, is solved numerically to obtain temperature 
profiles within the cell. Their model is used to show that as the rate of discharge 
is increased, the temperature increase is proportional to the increase in cell 
current. Szpak et al. (11) present a one-dimensional thermal model for disc- 
shaped Li/SOCl 2 cells, similar to the model by Parnell and Szpak (10), which 
accounts for catastrophic thermal runaway. This model is formulated by associating 
thermal runaway with ignition and burning, via one or more cell defects, rather 
than explosion. A differential energy balance and a differential material balance, 
which includes diffusion of species only, are solved simultaneously using a numerical 
method. The model is used to predict the time and position dependent temperature 
and concentration profiles. The model predictions compare favorably to what 
is physically observed which may indicate that the reactions causing thermal 
runaway are initiated by localized heat sources developed from defective cells or 
cell components. Cho and Halpert (12) have presented a simple thermal model 
for spirally wound Li/SOCl 2 primary cells. In this work, an energy balance which 
requires an experimentally determined term for the heat generation rate is solved 
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analytically to obtain the temperature of the cell as a function of time. They 
assume a uniform cell temperature which changes with time as has been done by a 
number of workers (8, 13-16). Cho and Halpert (18) and Cho (19) have taken the 
modeling one step further by developing a means to calculate the heat generation 
rate term in their model (12) using experimentally determined resistances. An 
electrical circuit analogue is used which consists of the thermal resistances of the 
various cell components. None of these models can be used here because each model 
treats the cell interior as one, pseudohomogenous material. 

Other Li/S0Cl2 cell models have been developed (13-17). Evans et al. (13) 
and Tsaur and Pollard (14-17) use conservation of mass and charge to determine 
the concentration and reaction rate profiles in the cell. Both use porous electrode 
theory (20, 21) and concentrated solution theory (22-24) to develop the governing 
equations in one dimension. Evans et al. ( 1 3) include four cell regions in their model 
whereas Tsaur and Pollard (14) essentially treat two cell regions. The boundary 
conditions used in each model are different and are based on different simplifying 
assumptions. It is assumed in each model that the cell temperature is uniform at 
any given moment in time. The cell temperature change with time is predicted 
using an overall energy balance. Tsaur and Pollard (15, 16) extend their model to 
include the additional species present when an acid electrolyte is utilized, such as 
in reserve cells (25). However, these previous models cannot be used to investigate 
the two-dimensional temperature profiles in spirally wound Li/SOCl 2 cells. 
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MODEL FORMULATION 


The differential energy balance for region k of a battery can be written as 


= V ■ (\ t VT) + Y,i"' 0 = -•?••••) (3] 

J 

where c Pi jt and A* are the heat capacity and thermal conductivity of region k , 
respectively. In some instances, the steady state tempera*. -are profiles are of interest 
in which case and the term dT/dt in Eq. [3] is set to zero. The q'j' terms in Eq. [3] 
represent various heat sources and sinks. Heat generated by polarization (q'”) and 
heat generated due to the entropy change of the current producing reactions ( q '"), 
reactions 1 and 2, will be included in Eq. [3] for normal discharge conditions. The 
heat generation terms q”' and q'j' are calculated as follows 26. 27) 


+ ?; 


gp + qa 

V ceI / 


[4] 


where 

q p + q a = i(e oc -E- = I(E tn - E) [5] 

E tn in Eq. [5], termed the thermoneutral potential, is the theoretical open circuit 
potential of the cell at absolute zero. The cell voltage. £ in Eq. [5], must be 
specified and can be obtained from experiment or predicted using an earlier model 
(13). Other sources/sinks of heat would include exothermic or endothermic chemical 
reactions occurring in the cell which do not produce current (27). These sources 
are assumed to be negligible in this paper as is usually done for normal operation. 
To simulate thermal runaway, an additional q'j' term is included in Eq. [3], qj', 
representing additional heating caused by, perhaps, a cell defect and/or exothermic 
chemical reactions. This heat source is assumed to be localized causing a hot spot 
in the cell. This approach to simulating thermal runaway was used by Szpak et 
al. (11). 
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To complete the mathematical description of the cell, the boundary conditions 
must be specified. For the cases investigated here, either a zero flux (insulated ' 
boundary condition 


, dT , dT 
Xk > z di Uz + Xk ' y dy Uy 


r 6' 


or a convection boundary condition 


. dT dT 

~ Xk ' z dx nz ~ Xk ’ y dy n 

is used. In Eq. [7], T w is the temperature of the cell at its outer wall and is a 
function of x, y, and t, Ta is the surrounding temperature and is a fixed value, 
and h is a heat transfer coefficient. It is assumed here, for simplicity, that h is a 
constant and adequately characterizes the heat transfer through the cell casing and 
from the case wall to the surroundings via convection. 


y = h(T w - T a ) 


J 
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SOLUTION PROCEDURE 


The finite element method (28, 29) is used to solve numerically Eq. [3], with 
Eq. [6] and [7], for the model region shown in Fig. 3. The computer package called 
TAPAZ2D (29) is used. Input files to T0PAZ2D, containing the coordinates of 
the finite element mesh, are created using another computer program called MAZE 
(30). 

In order to obtain reliable, consistent results using the finite element method, a 
smooth mesh must be constructed and nodal points at region interfaces must match 
up. To achieve this, the two-dimensional spiral mesh is developed by constructing 
two semicircular meshes, offsetting one semicircular mesh, and then joining the two 
meshes together. A schematic of a model region constructed in this way is shown 
in Fig. 3where x and y are the orthogonal vectors defining the two-dimensional 
plane and which intersect at the origin. The vector x x is shown in Fig. 3 because 
a temperature profile along this vector is shown later in this paper. The left half 
consists of four concentric semicircular cells and the right half consists of three 
concentric semicircular cells. Comparison of Fig. 2 and Fig. 3 shows that this 
construction provides a reasonable representation of the spiral design. For normal 
discharge conditions, meshes were refined until three digit accuracy was obtained. 

Table I lists the cell specifications for the spirally wound design investigated 
here. A similar cell is currently being developed by Wilson-Greatbatch (31) and 
may be used for applications in space (32). The thermophysical properties of each 
cell region and other model inputs are listed in Table II. Note that the properties 
for the cathode region and separator region axe average properties calculated by 
assuming that all voids are filled with electrolyte. Here, the porous cathode was 
assumed to be 85 % porous and the separator 70 % porous. 
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RESULTS AND DISCUSSION 


Two simplified cases of the spiral design were investigated. The spiral was 
approximated using concentric circular regions, the top schematic in Fig. 4. and a 
core region, the bottom schematic in Fig. 4. The core region approach is based on 
the assumption that the cell interior can be represented by one homogeneous region. 
This core region possesses thermophvsical properties which are the average of the 
properties of all the cell components. The core region properties were calculated as 
follows 

- Si 6 jPj 
P ~ Lifi 

where 6; is the thickness of region i and pi is the value of a thermophysical property 
of region i as listed in Table II. The core properties are listed in Table II. These two 
cases simplify the calculations because the governing equations need be solved in one 
dimension only. In fact, the steady state equation set for the core approach is easily 
solved analytically. Also, these two cases represent two different extremes. The 
circles case should yield a maximum temperature drop from the center of the cell to 
the cell exterior because all heat is forced to flow radially through all cell regions, 
including the separator region which offers the most resistance to heat flow. The 
core case should yield a minimum temperature drop from the center of the cell to the 
exterior because the thermophysical properties of the core region closely resemble 
those of the most conductive regions, as shown in Table II. The spiral case, as shown 
in Fig. 3, should yield temperature profiles somewhere between these two extremes. 
Heat flow in the radial direction is inhibited by the separator; however, heat can flow 
out the spiral along the more conductive materials. In the results which follow, the 
core approximation, the concentric circular regions approximation, and the spiral 
case are referred to as core, circles, and spiral, respectively. 
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Comparison of the core, circles, and spiral models: 


Figure 5 shows predicted time-temperature profiles for the core, circles, and 
spiral cases for a 4 A discharge. The heat generation rate program, shown in Fig. 6, 
was calculated using Eq. [5] and the voltage transient shown in Fig. 6. The voltage 
transient was approximated from experimental results reported in the literature (2, 
6. 33). Voltage transients can also be predicted from earlier models (13. 14). It 
should be noted that there exists a wide range of data for the Li/SOCh cell reported 
in the literature due to differences in the discharge rates, electrode materials, and 
cell designs investigated (2-6, 25, 33-36). Notice in Fig. 5 that the predicted 
maximum and minimum temperature dependence on time for the spiral case falls 
between those of the core and circles cases, as expected. Also note that the core 
case shows an essentially uniform cell temperature, whereas the circles and spiral 
cases show about a 2 K temperature drop from the cell interior (T max ) to the cell 
wall (Tmin). The general shape of the curves are in qualitative agreement with 
experimental temperature-time curves (6), as shown in Fig. 9; an initial rise in 
temperature as the electrodes are first polarized, then a plateau is reached where 
heat generation is about equal to heat transfer away from the cell, then a sharp 
increase in cell temperature at the very end of the discharge. A comparison of 
model predictions presented here and experimental results (6) is made below. 

The two-dimensional temperature profile for the spiral cell at approximately 
1.5 hr of discharge (corresponding to the dot in Fig. 5) is shown in Fig. 7. The 
plotted points are the temperatures at each node point of the spiral mesh. The 
temperature is plotted over the two-dimensional plane defined by the orthogonal 
vectors x and y shown in Fig. 3. The maximum temperature is in the electrolyte 
region located near the center of the cell (see Fig. 3). As shown in Table II, 
the electrolyte possesses a low thermal conductivity and thus serves to retain heat 
relative to the other cell components. The temperature drops from this maximum 
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to the minimum temperature located at the outer comer of the spiral, point A in 
Fig. 3. The largest temperature drops are across the separator and appear as the 
void spaces between the various levels of the spiral. 

Figure 8 compares spatial temperature profiles as predicted using the core, 
circles, and spiral methods. These profiles correspond to approximately 1.5 hr into 
the discharge, as shown by the box and dot in Fig. 5. Two profiles are shown for the 
spiral case, one along vector y and one along vector x\ corresponding to the vectors 
shown in Fig. 3. In order to better compare the profiles, the temperature difference 
T -Train is plotted on the ordinate. Figure 8 show's that the core approximation does 
not adequately predict the temperature drop in the spirally wound cell. In this case, 
the predicted temperature drop is about 2 K according to the spiral model whereas 
the core approximation yields only about a 0.2 K drop. Fig. 8 also show's that the 
profiles for the circles approximation and spiral case are in close proximity. The 
circles approximation is useful because it predicts maximum temperature profiles, 
as shown in Fig. 5 and 8. and is easily solved in one-dimension using TOPAZ2D. 

Several items are noteworthy in the forgoing comparisons between the core, 
circles, and spiral cases. The spiral model region, shown in Fig. 3, has boundaries 
different from the core and circles model regions, shown in Fig. 4. The difference 
is the exposed end of the spired cell between points A and B in Fig. 3. For the core 
and circles cases, Eq. [7] is used at the outer boundary. For the spiral, Eq. [7] is 
used for the outer boundary (clockwise from point A to point B in Fig. 3) and Eq. 
[6] is used for the end of the spiral cell (counterclockwise from point A to point B 
in Fig. 3). This use of insulated boundary conditions at the exposed end represents 
a worst case; in the actual cell, shown in Fig. 1, the lithium electrode is extended 
and the rest of the spiral end borders electrolyte. For comparison purposes, Eq. [7] 
w'as used for all boundaries of the spiral. The relative temperature drop from the 
cell interior to the exterior remained essentially the same, however, the absolute 
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temperature profile dropped about 1 K. What can be gleaned from this and the 
forgoing comparisons is that the spiral design aides in cooling down the cell, as 
opposed to the circles case, because heat can be conducted out the spiral path as 
well as out the radial direction. This increase reduces the cell temperatures but 
does not substantially reduce the temperature drop from the center of the cell to 
the cell wall. 

Finally, the effect of different heat generation rate distributions on the temper- 
ature profiles was briefly investigated using the spiral model. Four cases were run: 
90 % of the heat generation contained within the lithium anode and 10 % within 
the cathode, 10 % within the anode and 90 % within the cathode, 90 % within the 
separator and 10 % within the cathode, and 50 % within the anode and 50 % within 
the cathode. Each case yielded about a 2 K temperature drop. 

Comparison of circles model with experimental results: 

The one- dimensional circles model was used to simulate the temperature 
profiles for an experimental test case reported by Abraham et al. (6) for a C-size 
Li/SOCh cell. The component dimensions, cell volume, and other test conditions 
given in their paper were used in the model; eight wraps were assumed. The heat 
generation rate program was calculated based on the voltage transients shown in 
Fig. 7 of their work. The value of the heat transfer coefficient was set to 15 
W/m-K, instead of 10 W/m-K as in Table II, in order to match predictions to the 
experimental data as closely as possible. 

The experimental temperature-time curves for the 4 A discharge of a catalyzed 
Li/SOCl 2 cell are shown in Fig. 9. The difference in cell temperature from the 
center of the cell to the can wall is shown to be as much as 13 K. Model predictions, 
not those shown in Fig. 9, show only about a 3 K to 4 K temperature drop. 
The model predictions given in Fig. 9 show about an 8 K temperature drop and 
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were obtained by making several changes to the model as described belo%v. Several 
factors could account for the differences between the model predictions and the 
experimental measurements. First, this data is for a cell having porous cathodes 
containing 5 weight percent dibenzotetraazaannulene complex of cobalt Co-TAA). 
These catalyzed cells exhibit discharge behavior different from uncat alyzed cells (6) 
and the heat generation rate, as predicted using Eq. [5] and the thermoneutral 
potential listed in Table II, used in the model may not be appropriate. Assuming 
that this difference has no effect on model predictions, several other factors may be 
considered. The heat generation within the cell may not be uniform; it is probable 
that heat generation is greatest near the center of the cell and diminishes in value 
when approaching the cell ex'erior. This can be explained by considering the current 
density; near the center of the spirally wound cell there exists less electrode area 
over which the radially flowing current is distributed than at radii further from the 
center of the cell. Thus the electrochemical reaction rate would be greatest at the 
center, evolving more heat at the center, and decrease from there to the cell exterior. 
Decreasing thermal conductivities could also explain these large temperature drops. 
As the cell is discharged, the SOCI 2 electrolyte is reduced at the cathode and is 
sometimes depleted to the point where cell gases fill some of the voids previously 
occupied by electrolyte (32). To approximate these possibilities two changes were 
made to the model. The heat generation rate was set to a maximum value in the 
first concentric cell (relative value of 4) and was decreased linearly to a minimum 
value in the eighth concentric cell (relative value 1). The densities, heat capacities, 
and thermal conductivities of the regions containing electrolyte were changed by 
assuming that at the end of the discharge only SO 2 remained in the void portions 
of these regions. These properties were changed linearly throughout the discharge 
from their starting values, reported in Table II, to their ending values, reported 
in Table III. The model predictions are superimposed on the experimental data 
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in Fig. 9. Agreement between predictions and experiment has been improved by 
these changes. Steeper heat generation rate profiles were tried, resulting in larger 
temperature drops, but the temperature-time relationships did not agree well with 
the experimental results in Fig. 9. It should be possible to improve the fit of the 
theoretical predictions to the experimental results by using a parameter estimation 
technique in conjunction with T0PAZ2D. 

These differences between experiment and model indicate several possibilities. 
Perhaps a more detailed program for the heat generation rate is needed. Other 
chemical or electrochemical reactions could be occurring, especially at these high 
rates of discharge, which cause greater heat generation rates; these reactions could 
be localized. The thermal runaway simulations presented next were run with 
some of these thoughts in mind. More experimental data, including temperature 
measurements throughout the cell from center to wall, axe needed to analyze 
thoroughly the temperature profiles in the spirally wound cells and to refine the 
modeling. 

Thermal runaway simulation: 

Explosive failures are often associated with localized temperatures above the 
melting point of Li (2, 3). Therefore, initiation of thermal runaway was simulated 
by assuming that a hot spot, located near the center of the cell, is formed via 
a localized exothermic chemical reaction (for example, between Li and sulfur or 
a sulfur compound) or cell defect. Two simulations were run with the hot spot 
located at two different points in the cell. For case one, the hot spot was located at 
the center of the cell in the separator region bordering the cathode and electrolyte. 
For case two, the hot spot was located in the center of the electrolyte in the middle 
of the spiral. Case two was also approximated using the one-dimensional circles 
model for comparison purposes. The hot spot areas were approximately 0.031 mm 2 
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and 0.029 mm 2 for case one and two, respectively. The finite element meshes were 
refined until results were accurate to within 1 K for case one and 3.5 K for case two. 

The conditions under which thermal runaway usually occur include high 
discharge rates and high temperatures (2, 3-5). Therefore, the inputs to the model 
included an environmental temperature of 390 K, a uniform initial temperature in 
the cell of 390 K, and a uniform heat generation rate of 1 x 10 5 W/cm 3 . The hot 
spot heat generation rate, q was set equal to 5 x 10 8 \Y/cm 3 . The simulations 
were run until the minimum temperature in the cell exceeded the melting point of 
Li (453.6 K). It was assumed that when temperatures exceed the melting point of 
Li, molten Li bridges the gap between anode and cathode, forming an internal short 
circuit. Upon forming the short, it was assumed that q’J' would become infinitely 
large, the temperature would soar, and in a matter of milliseconds the cell would 
vent or explode (11). 

Figure 10 shows the maximum and minimum temperature-time curves for 
case one, case two, and the circles approximation to case two. The maximum 
temperature increases much more rapidly in case 2 than in case 1. In case 2, the 
hot spot is isolated from the more conductive regions of the cell by the electrolyte. 
In case 1, the hot spot borders the cathode region and the heat is better conducted 
out the cell via the spiral path. Note that the spiral geometry still improves the 
heat transfer in case 2 as is indicated by the higher temperatures predicted using 
the circles model. 

Figure 11 shows a two-dimensional temperature profile for case two approx- 
imately 470 seconds after initiation of the hot spot, corresponding to the box in 
Fig. 10. A temperature spike is located in the electrolyte region where the hot 
spot is centered. Note that minimum temperatures exist around most of the cell 
perimeter except for a portion of the cell located near the end of the spiral (i.e., 
between points A and B in Fig. 3). This rise in temperature at the spiral end is 
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due to the insulated boundary conditions applied at this boundary. In the actual 
cell (Fig. 1) the end of the spiral borders electrolyte (or cell gases if the electrolyte 
has been used up) and one would expect heat flow to be inhibited at this juncture. 
This result indicates that it is important how the end of the spiral is treated in the 
actual design. If possible, the end of the spiral should be in close proximity to the 
case can and void space should be minimized. Researchers (3) have observed that 
the temperature distribution over the exterior wall of the cell varies, some regions 
at high temperatures relative to the rest of the exterior surface. Perhaps these high 
temperature area s correspond to the location of the spiral end in the cell interior. 

Many situations, both normal operation and thermal runaway conditions, can 
be simulated using the model. These simulations can lead to a better understanding 
of the thermal behavior of this design as has been achieved here. Indeed the power 
of the model lies in the ability to investigate various scenarios simply by changing 
model inputs. The modeling can also be used as an aid in the interpretation of 
experimental data, as has been done here. 
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CONCLUSIONS 


A two-dimensional thermal model of a spirally wound Li/SOCU cell has been 
developed. The model and two one- dimensional approximations have been used 
to understand better the thermal behavior of this battery. Comparison of model 
predictions with experimental data support two contentions. The heat generation 
rate seems to be greatest at the center of the cell and diminishes towards the exterior 
of the cell. Also, thermal conductivities seem to decrease as the cell discharges. 
This could be due to the increasing volume occupied by cell gases as the discharge 
proceeds. Two thermal runaway situations were investigated using the model based 
on the theory that localized hot spots initiate thermal runaway. The effectiveness 
of the spiral design in conducting heat out of the cell is dependent upon where these 
hot spots are located in the cell. If hot spots are in contact with those regions of 
the cell having relatively high thermal conductivities then the spiral design serves 
to improve significantly heat dissipation. 

The model could be modified to investigate other physical situations. For 
example, forced convection, radiation, and conduction (via a potting material, 
e.g.) at the cell boundaries could be included in the boundary conditions of the 
model. Also, the model could be applied to other spirally wound batteries simply 
by changing model inputs (e.g., thermal conductivities, thermoneutral potential). 
Improvements to the model might include extending the model to three dimensions 
using the computer code T0PAZ3D (37), including contact resistances between 
regions, and adding temperature dependent thermophysical properties (10). 
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LIST OF SYMBOLS 


c p jt specific heat of material k , J/Kg-K 

E cell voltage, V 

E oc open-circuit voltage of the cell, V 

E tn thermoneutral potential of the cell. Y 

F Faraday’s constant, 96487 C/mol of electrons 

h heat transfer coefficient, W/m 2 K 

/ total cell current, A 

n r normal vector in the x direction 

n y normal vector in the y direction 

Pk property of material k 

p average property of core region 

q'" heat generation rate due to source j. W/m 3 

q'p heat generated due to cell polarization. W/m 3 

q'” heat generated due to entropy effects. W/m 3 

t time, s 

T cell temperature, K 

T 4 ambient temperature, K 

\ cel/ cell volume, m 3 

x dimension defining two-dimensional plane, m 

y dimension defining two-dimensional plane, m 

Greek Symbols 

6k thickness of region A:, m 

e porosity or void volume fraction 

A k thermal conductivity of region k, W/n-K 

Pk density of material k, Kg/m 3 

Subscripts 

A ambient conditions 

k region k 

j source j (either p, s, or d) 

oc open circuit 

w conditions at the wall 

x x dimension 

y y dimension 
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Table I. Cell specifications for a spirally wound Li/S0Cl2 cell 


Cell specification 

Value 

References 

Thickness of cathode 

1.12 mm 

31, 38 

Porosity of cathode 

0.85 

10, 38 

Thickness of anode 

0.7 mm 

31. 32. 35 

Thickness of current collector 

0.14 mm 

31, 32, 38 

Thickness of separator 

0.14 mm 

31, 32, 38 

Porosity of separator 

0.7 

38 

Number of wraps 

« 3 1/2 

31, 32 

D-size cell volume 

« 4.94 x 10 4 mm 3 

31, 32 



Table II. Model inputs 


Thermophysical properties 


Material 

Density 

(Kg/m 3 ) 

Thermal 

conductivity 

(W/m-K) 

Specific 

heat 

(J 'Kg-K) 

References 

Lithium anode 

534.0 

71.1 

3490.0 

10, 11, 39 

Nickel current 
collector 

8900.0 

90.5 

444.0 

11, 19 

Carbon 

1950.0 

23.8 

712.0 

10, 11 

Glass matting 

220.0 

0.242 

963.0 

10, 11, 38 

Electrolyte (1.8M 

LiAlCLt in SOCl 2 ) 

1690.0 

0.143 

1000.0 

11, 19, 38 

Cathode region 

1730.0 

3.69 

957.0 

10, 11, 19, 38 

Separator region 

1250.0 

0.173 

989.0 

10, 11, 38 

Core region 

1755.0 

23.82 

1577.0 

see Eq. [8] 


Model parameters 


Parameter 

h = 10 W/m 2 -K 
E tn = 3.723 V 


References 
11, 13, 14 
26 



Table III. Thermophysical properties of porous regions assuming that 
the voids are filled with SO 2 at approximately 20 psig and 60 °C 


Material 

Density 

(Kg/m 3 ) 

Thermal 
conductivity 
(W/m-K) * 

Specific 

heat References 

(J/Kg-K) 

S0 2 

5.54 

0.00816 

654.0 39 

Cathode region 

298.0 

3.577 

662.0 

Separator region 

70.0 

0.0783 

746.3 


LIST OF FIGURE CAPTIONS 

Fig. 1 Schematic representation of a spirally wound battery. 

Fig. 2 Top view of a spirally wound battery. 

Fig. 3 The two-dimensional model region for the spiral design. 

Fig. 4 One-dimensional approximations to the spiral design. The top 

schematic is the concentric circles approximation and the bottom 
schematic is the core region approximation. 

Fig. 5 Temperature-time curves predicted using the core, circles, 
and spiral models. 

Fig. 6 Heat generation rate and cell voltage for a 4 A discharge. 

Fig. 7 Predicted temperature profile after 1.5 hr of discharge. 

Fig. 8 Radial temperature profiles as predicted using the core, circles, 
and spiral models. 

Fig. 9 Comparison of predicted and experimental temperature- time 
curves for a 4 A discharge of a C-size Li/SOCl 2 cell. 

Fig. 10 Temperature- time curves for two thermal runaway cases. 

Fig. 11 Temperature profile in spirally wound Li/SOCL cell for thermal 
runaway case two approximately 470 s after hot spot initiation. 
Internal temperatures exceed the melting point of Lithium. 
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